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The emerging field of optomechanics seeks to 
explore the interaction between nanomechanics 
and light (see [lj for a recent review). Rapid 
progress in laser cooling of nanomechanical os- 
cillators [2j [3] promises new fundamental tests of 
quantum mechanics [4|, while applications bene- 
fit from ultrasensitive detection of displacements, 
masses and forces [5H7J. Recently, the exciting 
concept of optomechanical crystals has been in- 
troduced [8-10J, where defects in photonic crys- 
tal structures are used to generate both local- 
ized optical and mechanical modes that interact 
with each other. For instance, this opens the 
prospect of integrated optomechanical circuits 
combining several functions on a single chip (see 
also [11, 12]). Here we start exploring the collec- 
tive dynamics of arrays consisting of many cou- 
pled optomechanical cells (Fig. flk,b). We show 
that such "optomechanical arrays" can display 
synchronization and that they can be described 
by a modified Kuramoto model that allows to ex- 
plain and predict most of the features that will 
be observable in future experiments. 

The crucial ingredient of any optomechanical system 
is an optical mode (OM) whose frequency shifts in re- 
sponse to a mechanical displacement: 5uj op t = — Gx. 
This coupling, vice versa, gives rise to a radiation force, 
F = hG\a\ 2 , where \a\ 2 is the number of photons circu- 
lating inside the laser-driven OM. For a laser red-detuned 



from the OM (A = cj La 



<^opt < 0), dynamical back- 



action effects induced by the finite photon decay time n~ x 
lead to cooling of the mechanical motion. In contrast, 
for blue detuning (A > 0), anti-damping results. Once 
this overcomes the internal mechanical friction, a Hopf 
bifurcation towards a regime of self-induced mechanical 
oscillations takes place (Fig. lib) [13-18 . While the me- 
chanical amplitude A is fixed, the oscillation phase ip is 
undetermined. Therefore, it is susceptible to external 
perturbations. In particular, as we will see, it may lock 
to external forces or to other optomechanical oscillators. 

Synchronization has first been discovered by Huygens 
and is now recognized as an important feature of col- 
lective nonequilibrium behavior in fields ranging from 
physics over chemistry to biology and neuroscience [19], 
with applications in signal processing and stabilization 
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FIG. 1. Optomechanical crystals may be used to build ar- 
rays with several localized optical and mechanical modes, (a) 
Potential setup fabricated as a periodically patterned, free- 
standing dielectric beam on a chip with laser drive via tapered 
fibre as in [HJ[9]. (b) Schematic array of mechanically coupled 
optomechanical cells, (c) For a single cell, at sufficient laser 
drive power, there is a Hopf bifurcation towards self-induced 
oscillations with an undetermined phase (p. 



of oscillations. A paradigmatic, widely studied model 
for synchronization was introduced by Kuramoto [20] . 
For two oscillators, his phase evolution equation reads 
(pi = Vt\ + K sin((/?2 — <£i), and likewise (£2- One 
finds synchronization (ipi = £2) if the coupling K ex- 
ceeds the threshold K c = |^2 — fii|/2, and the phase 
lag Sep = (f2 — <Pi vanishes for large K according to 
sin(<5<£>) = (CI2 — Cti)/2K. For the globally coupled, mean- 
field type version of infinitely many oscillators, there is 
a phase transition towards synchronization beyond some 
threshold K c that depends on the frequency distribution 
[2T] , In many examples the Kuramoto model is found 
as a generic, reduced description of the phase dynamics. 
Nevertheless, for any specific system, it remains to be 
seen whether this model (or possibly a structurally sim- 
ilar variant thereof) applies at all, and how the coupling 
K is connected to microscopic parameters [221124] . We 
now turn to this question in the case of optomechanical 
oscillators. 

A single optomechanical cell consists of a mechani- 
cal mode (displacement x) coupled to a laser-driven OM 
(light amplitude a): 
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Here Q is the mechanical frequency, T the intrinsic damp- 
ing, G the optomechanical frequency pull per displace- 



ment, and a max is the maximum light-field amplitude 
achieved at resonance (set by the laser drive). 

Near the Hopf bifurcation (Fig. fit), we can capture 
the essential dynamics by eliminating the light field [16] 
and switching to the phase- and amplitude-dynamics of 
the resulting Hopf oscillator: 
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Here A is the steady- state amplitude, and 7 is the rate at 
which perturbations will relax back to A. Both depend 
on the microscopic optomechanical parameters, such as 
laser detuning and laser drive power, and both vanish at 
the bifurcation threshold (see methods section). More- 
over, we have introduced an external force F(t) (as added 
to Eq. ([2|) 

We start our discussion by considering phase-locking 
to an external force F(t) = Fq sin(ujpt). To this end we 
time-average Eq. ([3|, keeping only the slow dynamics, 
under the assumption ujp « Q. This results in 



Sep = —SQ + Kp sin((5(/?) . 



(5) 



where Sep = (p(t) + uopt, SO, = Q — ujp, and Kp = 
Fo/2mflA. Eq. ^ is a special case of the Kuramoto 
equation. Direct numerical simulation confirms the good 
agreement between the microscopic optomechanical dy- 
namics and the simplified descriptions via Eqs. (3J4) and 
([5|. In Fig. [2] we show sin 5 cp(t) and its time-average 
(sin 5 ip(t)), with the phase cp being extracted from the 
complex amplitude of motion, f3 = x + ix/Cl. Phase- 
locking sets in when Sep = has a solution, i.e. for 
\SCt\ < Kp, resulting in an "Arnold tongue" (see Fig.^t 



We now turn to the dynamics of coupled cells, each of 
which is described by Eqs. (HI [2]). To these equations, we 
add a mechanical coupling, set by a spring constant k: 
m'x\ = . . . + k(x2 — xi). In the Hopf model, this yields 
a force F\ = kA2 cos((/?2) on the first oscillator (and vice 
versa). The case of optical coupling will be mentioned 
further below. 

In order to arrive at the time-averaged dynamics for 
the phase difference, Sep = ip2 — <Pi, it is necessary to 
go further than before, keeping A(t) = A + SA(t) in the 
phase equation, and eliminating the amplitude dynamics 
to lowest order (see methods for the derivation; and [25] 
[26] for further examples where the amplitude dynamics is 
crucial). Then, we arrive at an effective Kuramoto- type 
model for coupled optomechanical Hopf oscillators: 



Sep = —SQ — C cos(Sip) — K sin(26cp). 



(6) 



In contrast to the standard Kuramoto model, 25cp ap- 
pears, which will lead to both in-phase and anti-phase 
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FIG. 2. Phase- locking of an optomechanical cell to an exter- 
nal force. (a,b) Phase lag S(p between oscillations and force, 
displayed via sin Sip, outside (a) and inside (b) the phase- 
locked regime (colored in (c)). Time is plotted in units of 
2tt/5£1. (c) Plot of (sin 5cp) as a function of frequency dif- 
ference 5Q, at the force strength indicated by the blue tri- 
angle in (d), compari ng op tomechanics (blue, solid) against 
the Hopf model (Eqs. ( |3|4| ); magenta, dash-dot) and the Ku- 
ramoto model (Eq. ([5]); black, dash), (d) (sin Sip) in the force 
vs. frequency difference plane, showing the Arnold tongue. 
Deviations between the different models set in at high forces. 
(Lines indicate the transition towards phase- locking, styles as 
in (c)). (Microscopic cell parameters: A/Q = 1, k/Q — 1, 
r/Q = 0.01, hG 2 a 2 max /mQ 3 = 0.36). 



synchronization. This corresponds to two distinct min- 
ima in the effective potential that can be used to rewrite 
Eq. ^: Sip = -U'(5tp). The coupling K = k 2 /2m 2 n 2 j 
diverges near the bifurcation, where 7—^0. In the fol- 
lowing we focus on the case of nearly identical cells where 
the coupling C can be neglected; C/SQ = k/2mVf? <C 1. 

To test whether these features are observed in the full 
optomechanical system, we directly simulate the motion 
and increase the coupling k for a fixed frequency differ- 
ence SQ — O2 — fii- The results are displayed in Fig.[3ja- 
c). Beyond a threshold fc c , the frequencies and the phases 
lock, indicated by a kink in (sin Sip). As the coupling in- 
creases further, the phases are pulled towards each other 
even more, so \Sip\ decreases. Thus, coupled optome- 
chanical systems do indeed exhibit synchronization. As 
predicted, there is both synchronization towards Sip —> 
and Sip — >> 7r. 

The dependence of the threshold k c on the frequency 
difference SQ is shown in Fig. |3]i. The observed behav- 
ior k c ~ \Z r SQ at small Sft is correctly reproduced by the 
generalized Kuramoto model, Eq. ([6|. For Sft > 7 devi- 
ations occur via terms of higher order in 5Q/j, starting 
with — (SQ/j)Kcos(2Scp) in Eq. (|6|. These produce a 
linear slope k c ex SCI, see Fig. [3]i. 

In terms of experimental realization, optomechanical 
crystals [8] E] offer a novel promising way to build ar- 
rays of optomechanical cells. They are fabricated as 
free-standing photonic crystal beams (Fig. Ilk). Varia- 
tions of the /im-scale lattice spacing produce both local- 
ized optical and vibrational modes. The strong confine- 
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FIG. 3. Phase-locking for two coupled optomechanical cells, (a) The phase particle in the effective Kuramoto potential U(5ip), in 
the de-synchronized and phase-locked regime (left/right), (b) Time-evolution of sin (Sip (t)). (c) When the mechanical coupling 
k exceeds a threshold, the phase difference 5(p between the oscillations locks in spite of different bare mechanical frequencies, 
here SQ — 0.003 Oi. Both in-phase and anti-phase regimes are observed, (d) Phase lag, expressed via (sin(£<£>)), in the plane 
coupling k vs. frequency difference 5Q, including a comparison of the critical coupling k c (white, solid) with the one from a 
Hopf model with one fit parameter (green, dash) and the generalized Kuramoto model (blue, dash-dot). (Cell parameters as 
in Fig. [2}. 



merit leads to extremely large optomechanical couplings, 
on the order of G ~ 100 GHz/nm. Typical parameters, 
that we use in our simulations for experimentally realis- 
tic results, are G = 100 GHz/nm, mechanical frequency 
Q = 1 GHz, mass m = 100 fg, mechanical quality fac- 
tor Qm = Q/T = 100, cavity decay rate n = 1 GHz 
and laser input powers such that the circulating photon 
number |a max | 2 ^ 100. 

To consider optomechanical arrays like in Fig. [TJi, we 
use finite element methods (FEM) to simulate two iden- 
tical cells arranged on the same beam (Fig. ELb). The 
optical and vibrational couplings mediated by the beam 
can be deduced from the splitting between the resulting 
symmetric/antisymmetric modes. The results shown in 
Fig. [4J3, d validate the parameters considered above and 
indicate mechanical couplings k/mVt 2 < 0.01. Due to 
the relatively strong optical coupling (~ THz), distinct 
OMs in the individual cells can only be achieved by pat- 
terning them to have frequencies sufficiently different to 
prevent hybridization (Fig, [4^). This requires different 
laser colors to address each cell. 

In experiments, a convenient observable would be the 
RF frequency spectrum of the light intensity emanating 
from the cells, \a\ (lj). We first show the spectrum as a 
function of frequency difference SQ for two mechanically 
coupled cells, driven independently (Fig. [5k). Frequency 
locking is observed in an interval around SQ = 0. Exper- 
imentally, the mechanical frequencies can be tuned via 
the "optical spring effect" [8]. 

The most easily tunable parameter is the laser drive 
power (oc c^max)- Synchronization sets in right at the 
Hopf bifurcation. For two cells (Fig. |5fc), we recover 
the regimes of in-phase and anti-phase synchronization. 
They differ in the synchronization frequency, Q(tt) — 
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FIG. 4. FEM simulation for two coupled cells localized on the 
same beam in an optomechanical crystal setup (cf. Fig. fla). 
(a) Horizontal displacement of vibrational "pinch" modes, (b) 
Transverse electric fields of optical modes. (c,d) The inter- 
cell couplings decay exponentially as a function of distance. 
(e) Blue-detuned lasers could illuminate two distinct optical 
modes (top) or a single hybridized "molecular" mode. (Geo- 
metrical parameters see methods section). 



£l(0) = k/mfl. At higher drive, we find a transition 
towards de-synchronization. This remarkable behavior 
can be explained from our analytical results. We know 
that 7 increases away from the Hopf bifurcation (i.e. for 
higher drive), leading to a concomitant decrease in the 
effective Kuramoto coupling K oc I/7, and finally a loss 
of synchronization. Near the transition, the frequencies 
fan out (as 5Q e ff °£ V<^max — ol c ). The multitude of peaks 
is produced due to nonlinear mixing. In some regimes, 
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FIG. 5. Mechanical frequency spectra of optomechanical arrays. (a,b) Frequency locking upon changing the detuning SQ = Q2 — 
Qi between the mechanical frequencies (magenta, dashed) of two independently driven, mechanically coupled optomechanical 
cells (k/mQi = 0.015). (a) Spectrum I(uj) of intensity fluctuations, I(i) = |ai(£)| 2 + |c^2(t) | 2 , from an optomechanical model. 



(c) Spectrum vs. laser input power for two independently driven, mechanically 
(d) Example of trajectories Xi(t), in units of k/G, displaying strong amplitude 



(b) Spectrum of x\ + X2 in the Hopf model, 
coupled cells (k/mfll = 0.01, 6Q = 0.005Qi) 
modulation, not described by the Kuramoto model (at a power indicated by the blue triangle in (c)). (e) Spectrum I(uj) vs. 
laser input power for an array of five mechanical modes (k = 0), coupled to a common optical mode. At large drive, a regime of 
chaotic motion is entered. Due to the presence of multiple attractors, the regimes observed in such a diagram may depend on 
how the parameters are swept, (f) Example of a trajectory in the chaotic regime (blue triangle in (e)). (color scale in all plots 
indicates |/(cj)| in units of the pea k he ight at uj — for a system at rest; 5 peaks are broadened for clarity; cell parameters are 
Ai = Ki — 100r^ = f^i, as in Figs.|2|3[ yellow triangles indicate the bare mechanical frequencies Qi) 



the Kuramoto model fails (Fig. |5j: 

For large arrays, it will be most practical to have iden- 
tical OMs that combine into extended 'molecular' modes, 
one of which is then driven by a single laser via an evanes- 
cently coupled tapered fibre (Fig. [la, Fig. He). For effi- 
cient excitation of self-induced oscillations, one has to 
ensure that the detuning A is equal to all the mechani- 
cal frequencies Qj in the array, to within |A — Clj\ < n. 
For arrays of reasonable size, the splittings between ad- 
jacent optical molecular modes will be more than 10 2 
times larger than k, such that we can ignore all but one 
OM. This setup then leads to a global coupling of many 
nanoresonators to a single OM, such that 



a 



i(A + E G i 



(7) 



and the force on each resonator is given by — HGj\a\ 2 . 
This setup comes close to realizing the all-to-all coupling 
most often investigated in the literature on the Kuramoto 
model. 

For illustration, we chose N = 5 cells (Fig. |5fe). As 
before, we find synchronization regimes. In addition, at 



higher drive, a transition into chaos takes place. Analyz- 
ing the time-evolution in more detail, we observe tran- 
sient fluctuations in amplitude and phase, with a strong 
sensitivity on changes in initial conditions (Fig.lEt). Note 
that in a single optomechanical cell one may also find 
chaotic behavior [14], but for far larger driving strengths. 



Optomechanical arrays open up a new domain to study 
collective oscillator dynamics, with room-temperature 
operation in integrated nanofabricated circuits and with 
novel possibilities for readout and control, complement- 
ing existing research on Josephson arrays [22 , laser ar- 
rays [23] and other nanomechanical structures [24] [27] . 
Recent experiments on 2D optomechanical crystals [28] 
could form the basis for investigating collective dynam- 
ics in 2D settings with various coupling schemes. Ap- 
plications in signal processing may benefit from phase 
noise suppression via synchronization [29]. Variations of 
the optomechanical arrays investigated here may also be 
realized in other designs based on existing setups, like 
multiple membranes in an optical cavity [30] or arrays of 
toroidal microcavities [2j [5] . 
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Methods 

The dynamics of a single optomechanical cell, 
Eqs. (TTl [2]) , in the self-induced oscillation regime can 
be mapped to a Hopf model close to the Hopf bifurca- 
tion. This model is described by the steady state ampli- 
tude A and amplitude decay rate 7 (see Eqs. (|3j [I])). 
The dependence of 7, A on the microscopic parame- 
ters can be deduced by expanding the average mechan- 
ical power input provided by the radiation pressure 
force, hG(\a\ x) (see Eq. (|l|), in terms of a = GA; 
HG(\a\ 2 x) = (fto4 ax 7r 2 )a 2 + (ha^ir^/O 2 )^ + 0(a 6 ). 
This yields 



and 



7 = 2TO7T2 - r 



(GA/nf =7^/(-27T 4 ^ 2 ), 



where the dimensionless coefficients 7r 2 (A/fi, k/0) and 
7T4(A/f2, k/O) only depend on the rescaled detuning and 
cavity decay rate. V = hG 2 a 2 nax /m0 3 is the rescaled 
laser input power. Numerically, the Hopf parameters 
may be found (even away from threshold) by simulating 
the exponential relaxation of the cell's oscillation ampli- 
tude towards A, after a small instantaneous perturbation 
of the steady-state dynamics. In general, to compare op- 
tomechanics to results from Hopf (e.g. Fig.^j, the optical 
spring effect also needs to be considered. 

Two mechanically coupled optomechanical cells are 
modeled as distinct Hopf oscillators with phase dynam- 
ics (pi(i), if2(t) and amplitude dynamics Ai(i), A 2 (t) 
according to Eqs. pi) and Q. The coupling forces 
read F\ = kA2 cos((/? 2 ), F 2 — kA\ cos(y?i). With 
Ai(t) = Ai + 5Ai(t), the solution for the amplitude dy- 



namics is SAi(t) = J_ e 



-7(t-t') 



fi(t f )dt' where h (t) = 



m Q 2 COS( P2(t)sm(fi(t) and likewise for f2(t). In 
the following, we consider small couplings, where SAi <C 
Ai. Then 5A\ (and likewise SA 2 ) is found to be 

5A ± (t) k 



A x 



2m\0\ 

e i(<P2(t) + <p!(t)) 

Im ' 



p i((p2(t)-(p!(t)) 



7-i(0 2 H-0i) 7-i(0 2 -Oi) 



Thus we can eliminate the amplitude dynamics to 
lowest order from the phase equations, by expanding 



A 2 (t)/A 1 (t) « A 2 /A 1 + SA 2 /A 1 - A 2 SA 1 /A 1 in the fol- 
lowing equation (likewise for (p 2 ): 



V\ 



-o x 



k 



miQi 



A 1 



cos((/? 2 )cos((pi). 



We now perform a time average, keeping only the 
slow dynamics near frequencies and ±\0 2 — ^i|. 
This leads to the stated result for the effective Ku- 
ramoto model, Eq. (|6]), after setting Sip = (f 2 — 
if\. The coupling constants are given by C = 
(fc/2) (A 2 /m 1 n 1 A 1 - A 1 /m 2 n 2 A 2 ) and K = (1_+ 
(6/6 +6/6)/2)A: 2 /4mi^im 2 ^ 2 , where Q = m^jA). 

The optomechanical simulation in Fig. [3] shows results 
for experimentally realistic microscopic parameters us- 
ing an input power well above the bifurcation threshold. 
This allows to observe the essential features predicted 
from Hopf and the effective Kuramoto-type model in an 
appropriate range of frequency detuning 50, . However, 
to achieve quantitative agreement of the Hopf model in 
Fig.[3j its parameter 7 has to be treated as an adjustable 
parameter (here 7 = 0.02^). Each simulation initially 
starts with a system at rest and considers an instanta- 
neous switch-on of the laser input power. Whether the 
system synchronizes towards Sep — >• or Sep — >• it also 
depends on the initial conditions. 

In principle, an amplitude dependence of the mechan- 
ical frequency, 0(A) ~ 0(A) + (dO/dA)SA can yield ad- 
ditional terms and an alternative synchronization mech- 
anism for two Hopf oscillators. However, numerical sim- 
ulations verified that for optomechanical cells this aspect 
can be neglected. 

For the finite-element simulation in Fig. [4] the unit cell 
in the periodic part is a l,396nm-wide, 362nm-long rect- 
angle with a co-centric rectangular hole of 992nm width 
and 190nm length. The thickness of the beam is 220nm. 
The isotropic Young's modulus of 168.5 GPa and the re- 
fractive index is 3.493. Each defect is 15 units in length, 
symmetric across the 8th(central) cell. The lattice con- 
stants vary linearly from 362nm at the edge to 307. 7nm 
at the center. The holes in the defects stay co-centric 
with the unit cell and remain constant in size. For a sin- 
gle cell these parameters have been reported in Ref. [9] . 
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